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An independent pair ansatz is developed for the many body wavefunction of dilute Bose systems. 
The pair correlation is optimized by minimizing the expectation value of the full hamiltonian (rather 
than the truncated Bogoliubov one) providing a rigorous energy upper bound. In contrast with the 
Jastrow model, hypernetted chain theory provides closed-form exactly solvable equations for the 
optimized pair correlation. The model involves both condensate and coherent pairing with number 
conservation and kinetic energy sum rules satisfied exactly and the compressibility sum rule obeyed 
at low density. We compute, for bulk boson matter at a given density and zero temperature, (i) the 
two-body distribution function, (ii) the energy per particle, (iii) the sound velocity, (iv) the chemical 
potential, (v) the momentum distribution and its condensate fraction and (vi) the pairing function, 
which quantifies the ODLRO resulting from the structural properties of the two-particle density 
matrix. The connections with the low-density expansion and Bogoliubov theory are analyzed at 
different density values, including the density and scattering length regime of interest of trapped- 
atoms Bose-Einstein condensates. Comparison with the available Diffusion Monte Carlo results is 
also made. 
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I. INTRODUCTION 

Models for interacting bosons, a focus of interest for 
decades, have gained renewed importance due to the 
intense current activity in Bose-Einstein condensates 

§11!- 

The Bogoliubov model for weakly interacting bosons 
H ^, Q has long been the starting point for more de- 
tailed theories [g, ||, [l^, [ll], [l2|. The ground state is 
approximated by an ideal Bose gas, and number conser- 
vation is violated by anomalous operator averages of the 
zero momentum condensate operator. Scattering out of 
(and depleting) the zero-momentum ground state pro- 
duces opposite-momentum pairs that are treated as ex- 
citations with a linear spectrum pj. 

Mean field pairing theories of superfluidity, by analogy 
with superconductivity, considered BCS-like anomalous 
pair-operator averages [Hj, preventing an undesir- 
able BCS gap in the excitation spectrum by a condition 
on the chemical potential . When neutron scattering 
estimated that in He II, the condensate fraction was less 
than 10% ||, [l| while the superfluid fraction is 100%, 
such models were given added impetus, with pairing as 
a part (or whole) of the ground state. More detailed 
many-body ground state calculations based on Jastrow 
wavefunctions go beyond the simple mean-field pairing, 



but to handle hard core potentials, must perforce have 
not only pairing, but also triplet, quadruplet and higher 
correlations. Thus many-body schemes like the hyper- 
netted chain (HNC) approximation become open-ended 
hierarchies of coupled equations @ . 

On the other hand, the discovery of Bose-Einstein con- 
densate (BEC) with a finite number of atoms prompted 
investigation of number-conserving reformulations of the 
Bogoliubov model, without invoking anomalous averages, 
in modified Bogoliubov excitation operators [Q, |l6). Al- 
though the BEC depletion fraction is estimated to be 
small, it is nonzero, and issues remain, of a generalized 
understanding of the 100% superfluid condensate Q. 

Motivated by low carrier-density superconductivity in 
high T c materials, negative-?/ fermion pairing models 
have shown (lTI Jl8| that real space pairing (as in pre- 
BCS models j9[, |l9[) occurs at low densities, and goes 
over to momentum space BCS pairing at high densities. 

In this paper, we consider a non- Jastrow ground state 
of number-conserving pairing that ignores higher correla- 
tions, but includes real space pairing correlations exactly. 
What we lose in a description of normal-state correla- 
tions and densities close to crystallization, we gain in a 
focus on the correlations responsible for superfluidity at 
lower densities. The emergent picture is of superfluidity 
from condensate plus real space pairs, with the relative 
fraction dependent on density and interaction, but the 
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superfluid fraction always unity at T = 0. Within a 
pair-only model, the hypernetted chain scheme for the 
ground state is now a closed form exactly soluble set of 
equations, rather than being a correlation hierarchy, to 
determine the two-body distribution function and the 
density matrices. Numerical solution of the Euler equa- 
tion, resulting from the Ritz variational principle, yields 
equation of states for any given soft-core potential, and a 
momentum distribution with the correct long wavelength 
behavior, having the proper Gavoret-Nozieres singularity 

Correlation functions show off diagonal long-range or- 
der (ODLRO), factoring into pairing functions that de- 
pend on the pair correlation defining the independent- 
pair ground state. Number conservation and other 
important sum rules are satisfied ^ The 
Bogoliubov-type results are recovered not at low den- 
sities, but in the limit of large pa 3 , where a is the inter- 
particle scattering length. 

The interaction-induced presence of nonzero- 
momentum particles yields an occupation number 
dependent wavevector that has some analogies with 
the ideal Fermi gas: one can define a nonzero Base 
wavevector fcf,, an average Bose energy Ei,, and a 
corresponding Bose temperature that all vanish as a 
function of scaled density/interaction, i.e. in the ideal 
Bose gas limit. 

More in detail, in this paper we use variational theory 
and Hyper Netted Chain methods ]?], ^2|, to add pair 
correlations to the mean field wave function, which is at 
the basis of the Gross-Pitaevsky approach. We consider 
the following ground state trial function which includes 
independent pair correlations (IPC) only 

*iP C (ri,...,r N ) = 1+J2 h (nj) 

+ 5? HnjMrim) + ..(;) 

' (i<i)^G<m) 

where the indices of the correlation factors h(rij), ap- 
pearing in the above summmed products, at all orders, 
never overlap. Writing h(r) in Fourier space, one can eas- 
ily verify that with the IPC trial function of Eq. (|l|) the 
zero-momentum condensate fraction interconverts with 
depletion pairs with zero total momentum only. Note 
that h is real i.e. condensation phase is locked to coher- 
ent pairing phases. Thus the ^ipc is characterized by a 
condensate and a coherent depletion. 

Although this condensate plus pair (only) ansatz is 
in the broad spirit of the Bogoliubov model there are 
three crucial differences: i) Number conservation is main- 
tained, with no anomalous averages; ii) The condensate 
and pairs are coherent and both are part of the ground 
state rather than the latter being excitations; iii) The 
expectation valule of the total hamiltonian is evaluated, 
without truncation, so the energy is a rigorous upper 
bound, and the correlation h(r) obtained by minimiza- 



tion goes beyond the Bogoliubov approximation. Exci- 
tations are not considered here, but will be related to 
density variations of the interaction-induced state of co- 
herent condensate-pair interconversion. 

The IPC model constitutes the underlying wave func- 
tion of Coupled Cluster theory, in the SUB(2) approxima- 
tion, which has been succesfully used in a number of stud- 
ies of strongly correlated systems j24| [2f| . In compari- 
son with the Jastrow ansatz, ^j(R) — IL<j(l + M r ij)); 
the correlation products h(rij)h(ri m ) ... in Eq. ([!]) with 
7^ (lm) constitute only a subset of those produced 
by tyj(R), and, as a result, ^>ipc(R) cannot be put 
in a product form like ^fj(R). Nonetheless, ^ipc(R), 
still satisfies the separability condition, namely for any 
{R} = {Ri,R 2 } Vipc(Ri,R2) -> ^ipc(Ri)^ipc(R2), 
in the limit of the subset of particles {Ri} being macro- 
scopically far from {i?2}- 

The trial function ^ipc(R) of Eq. (|l|), not being of the 
product form, cannot deal directly with two-body poten- 
tials which have strong singularities at short distances, 
like the hard sphere potentials or those of the Lennard- 
Jones type. However, the detailed shape of the two- 
body potential may not be relevant for low-density sys- 
tems, and suitable two-body pseudopotentials fl26j can 
be found to study in a quantitative way the dynamical 
and superfluid properties of Bose systems. 

We compute the two-body distribution function and 
one- and two-particle density matrices of a bulk boson 
system by using a new HNC cluster technique, based on 
the Renormalized Fermi Hyper Netted Chain (RFHNC) 
theory jp3|. This new cluster technique is made neces- 
sary because of the particular form of ^ipc(R), which 
requires the summation of reducible and unlinked dia- 
grams, contrary to the case of ordinary Jastrow theory. 
This novel method is based on treating the boson system 
as a collection of Fermi particles, having the same mass 
and interaction of the true bosonic ones, and a flavour 
degeneracy which equals the number of particles, so that 
they can all occupy the lowest state k — 0. The above 
method relies on the property of the FHNC cluster ex- 
pansion to be exact at any order in 1/N |2^, ^3|. This 
could later permit the study of finite-atom trapped con- 
densates. 

The plan of this paper is as follows. In Section II, we 
present the diagrammatics of the hypernetted chain ap- 
proximation, and show how the independent-pair model 
for the ground state allows for a closed set of equations. 
In Section III we relate the correlation functions (reduced 
or traced density matrix elements) to the pair correlation 
h(r) entering the ground state. In Section IV we derive 
the Euler equations for the optimal h(r) and the corre- 
sponding momentum distribution n(q). Results for dif- 
ferent soft-core potentials, and plots of various physical 
quantities are given in Section V. Finally, Section VI is 
a summary and conclusions. 
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II. HNC THEORY FOR THE INDEPENDENT 
PAIR MODEL 

We present in this Section the calculation of the two- 
body distribution function of a bulk boson system de- 
scribed by the IPC trial function of Eq. (|]). 

The pair distribution function 52(1,2) is one of the 
fundamental quantities in the study of interacting many- 
body systems. Both the structure function S(k) and the 
expectation value of the hamiltonian E[h] are given in 
terms of it. Therefore, its calculation is a necessary step 
for any quantitative analysis of the static and dynamical 
property of an interacting quantum many-body system. 
It is defined as 



i j 




ff2 (l,2) = g(r 12 ) 

N{N - 1) / dr 3 dr 4 . . . dr N \* IPC (R)\ 2 
P 2 fdR\y IPC (R)\ 2 ,( ' 

where the first equality is due to the translational in- 
variance of the system. The static structure function 
S(k) is strictly related to the pair function and given by 
I 

S(k) = l + pj dre*~ k %(r)-l) . (3) 

The two-body distribution function g(r) satisfies the 
following normalization sum rule M 



5(0) = l + p dr{g{r) - 1) = . (4) 



A. The new HNC scheme for the IPC model. 

A new many-body method, based on the properties 
of the FHNC cluster expansion §| and on the RFHNC 
summation technique p3| , is developed and used in the 
calculation of g{r\2). We show below, as done in Ref. 
H for the Jastrow model, that RFHNC applied to N- 
fiavour fermions (which are equivalent to a system of N 
bosons) produces the exact counting coefficients of the 
cluster diagrams in the independent pair ansatz. 

The standard HNC cluster decomposition is based on 
(i) a series expansion of the r.h.s. of Eq. (||) around 
h(r) — and ( ii) a diagrammatic representation of the 
various terms Kjfl , The independent pair condition in 
^ipc(R) implies that each vertex of any cluster diagram 
has at most two correlation lines reaching it, one com- 
ing from the ket l^/pc) (K-correlation) and the other 
coming from the bra (vP/pcl (B-correlation) . This prop- 
erty, on the one hand, kills all the Jastrow elementary 
diagrams, but, on the other, leads to a highly non trivial 
cluster expansion with extra vertex corrections that can 
be summed. A new pairing-focussed HNC treatment, 



FIG. 1: Diagrams showing the cancellation process of the 
HNC and the FHNC cluster expansions as disussed in the 
text. All the four diagrams give the same cluster integral 
I, and are constituted by the two substructures T(l, . . . , s), 
with s — 8, and 7(1, j). Diagrams (a) and (b) are bosonic. 
Diagrams (c) and (d) are fermionic, and their oriented lines 
denote the exchange function l(kpr). 

different from the standard one developed for Jastrow 
theory is needed to sum up the resulting cluster series, 
and (unlike Jastrow) it can be written in closed form. 

A first important subdivision of the cluster diagrams 
is between linked and unlinked diagrams, where the un- 
linked diagrams are those with two or more non overlap- 
ping parts. The linked diagrams are further subdivided 
in irreducible and reducible diagrams. Reducible (or sep- 
arable) diagrams correspond to cluster integrals which 
can be factorized in two or more parts. 

In standard HNC, all the unlinked diagrams cancel 
each other, as schematically shown by the following sim- 
ple example. The upper part of Fig. [j] displays two 
cluster diagrams, (la) and (lb) from different terms 
of the Mayer expansion. Both arc constituted with 
two linked diagrammatical structures, r(l,...,s) and 
7(1, j) = h(rij), where i or j may or may not coincide 
with one of the labels of T. The structure T(l, . . . , s) has 
s vertices, two of which correspond to the two external 
labels 1 and 2 of g{r\2)\ j) has two vertices only. Dia- 
gram (la) is unlinked, whereas diagram (lb) is linked but 
reducible. Therefore, their corresponding cluster terms, 
in Figure 1, have the same integral form 

I = J df x ... df s T(fi, ...,f s ) J dfidfjj(ri,fj). (5) 

The unlinked diagram (la) comes from the unlinked 
term ([r + 7]) in the numerator of Eq. (||) or as a product 
([Fx 7]) of the term (r) in the numerator and the term (7) 
in the denominator. The reducible diagram (lb) comes 
from a numerator term. Let us now compute for each 
of them the counting coefficient due to the presence of 
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j(i,j), under the assumption that the diagrammatical 
structure 7 can be attached to any of the s vertices of T 
in diagram (lb) (as in the case of Jastrow theory) EG] 



[r + 7l 
[r x 7 ] 



(N — s)(N 



1) 



2n 2 

N(N - 1) 
s(N - s) 



(6) 



Summing up the various coefficients, one can see that 
the leading terms (those in N 2 /il 2 ) cancel out. More- 
over, the next to leading order terms (those in N/fl 2 ) 
also vanish. This means that the combined diagrammatic 
structures (la) and (lb) do not contribute to the two- 
body distribution function in the thermodynamic limit. 
We have recovered in this example the well known prop- 
erty of the Mayer cluster expansion, that only irreducible 
diagrams contribute to g{r), and the reducible and un- 
linked diagrams correspond to terms of order 1/N or 
lower, which therefore vanish in the thermodynamic limit 

In the case of the IPC ansatz, not all the s reducible 
diagrams of type (lb) are allowed, as assumed in the pre- 
vious equation. For instance, if a given vertex, amongst 
the s possible ones, has both a B and a K correlation 
reaching it in the Y structure the second diagrammatical 
structure 7 cannot be attached to it (so e.g. lb can not 
occur). In this case, we do not have anymore the can- 
cellation of the next to leading order terms and one has 
to include I of Eq. (|J) (and, therefore both the unlinked 
and the reducible diagrams (la) and (lb)) in the cluster 
expansion with the proper coefficient. 

Let us now suppose that no reducible diagrams of the 
type (lb) are allowed. In this case, the leading coefficient 
of /, resulting from the sum of the contributions [r + 
7] a and [r x 7] a in Eq. (||), is proportional to N and 
given by —ap/Cl. The translational invariance of 7 brings 
in a factor 0, and the internal counting coefficient of T 
another factor tr x N(N - 1) . . JN - s + l)/il s , where 
tr is the topological factor of T pq] . Therefore the total 
coefficient, /(bosonic), of the cluster term / is non-zero 
and given by 



/(bosonic) = —tr X sp s 



(7) 



Unfortunately, the summation of reducible diagrams 
is very difficult in the general case, because it requires a 
direct counting of terms coming from the numerator and 
the denominator. 

The situation is much simpler in the Fermi case, where 
one can use the FHNC cluster expansion. With fermions 
of v flavours, because of the Pauli principle, the various 
expansion terms in both the numerator and the denomi- 
nator have coefficients, simply given by t x (— ) s+ 'p s /^ s ~', 
where t is a generic topological factor, s the number of 



vertices, I the number of Pauli loops (counting also the 
one-vertex loops). As an example, the terms with no 
exchanges have only one- vertex loops, therefore s = I, 
and the coefficient is t x p s . The vertices in a Pauli loop 
are connected by the two-body Fermi function l{kpTij), 
given by 



l(x) 



:(sin(x) — xcos(x)) 



(8) 



where x = kpr^ and kp — (6tt 2 p/is) 1 / 3 is the Fermi 
momentum. 

One may easily verify that, given the structure of the 
coefficients of the cluster diagrams, there is a complete 
cancellation of the denominator with the unlinked terms 
of the numerator |22[ . Therefore, referring to the example 
of Fig. [l], the corresponding fermionic diagrams of the 
cluster term I cannot be unlinked, as diagram (la). The 
fermionic diagram corresponding to (lb) is diagram (lc), 
where Pauli one-vertex loops have been added up. The 
loops, however, do not change / because l(kpru) = 1. Its 
global counting coefficient is tr X p s+1 . 

In addition, there are other reducible fermionic dia- 
grams, like diagram (Id), which have the same integral 
form of /. This is due to the normalization and com- 
pleteness properties of the Fermi function l(kprij) 



- / dfl(kpr) = - I drl 2 (k F r) = 1 , 



dr 3 l(kFri2)l(k F r 3 2) = l(k F r 12 ) 



(9) 



It is well known that for the Jastrow-Slater model 
all the reducible diagrams of the two-body distribution 
function cancel out exactly, as for the Bose case |^3| . One 
may understand such cancellation by looking at the two 
cluster structures (lc) and (Id). The coefficients of the 
two corresponding cluster terms are tr X sp s+1 for dia- 
grams of type (lc) and — tp x sp s+1 jv for diagrams of 
type (Id) (the separability point in both the structures 
(lc) and (Id) can be any of the s vertices of T). From 
the properties of the exchange function /(fcj?r) given in 
eqs. ([)]), it follows that the two cluster terms cancel each 
other. 

However for the IPC model, diagrams (lc) are not al- 
lowed for the same reason that (lb) is forbidden in the 
boson case (no B,K and 7 lines from a single vertex). 
Hence one is left in the IPC/FHNC case with uncan- 
celled diagrams of the (Id) type, and usingloop integrals 
like [I ends up in the cluster term / of Eq. (||) with a total 
coefficient 



/(fermionic) = — tr x sp 



s+1 



(10) 



But this is exactly as in the case of bosonic diagrams 
(la) and (lb): the FHNC result is the same as the HNC 
for bosons. 
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From the above discussion it follows that a conve- 
nient way to compute the two-body distribution function 
.9(712) of Eq. (||) is to consider a notional Fermi system, 
in which the fermions have the same mass of the origi- 
nal bosons, interact with the same potential, and have a 
flavour degenaracy equal to the number N of particles, 
so that they all occupy the lowest state with k — 0. We 
can then use the much simpler FHNC theory, under the 
constraint v = N, in a bosonic limit. 

This v = N limit implies l(kprij)/v — > 1/N. As a con- 
sequence, all the irreducible fermionic diagrams will van- 
ish in the thermodynamic limit, except those having one 
vertex loop only, which coincide with the corresponding 
bosonic diagrams. The reducible fermionic diagrams sur- 
vive, because integrations on the uncorrelated fermionic 
bonds bring Q factors. In the following, we will show how 
to calculate them in a simple and efficient way. 



B. Vertex corrections. 

In order to include reducible diagrams we use RFHNC 
theory which treats irreducible renormalized diagrams 
characterized by vertex corrections for each of their ver- 
tices ||. 

There are two kinds of vertex corrections 

• bosonic vertex correction Vf,. This vertex correction 
applies to vertices of type either B or K, which are 
reached only by one B- or if-correlation in the 
underlying irreducible structure; 

• fermionic vertex correction Vf. This vertex cor- 
rection applies to vertices of the type KB, namely 
those reached by both K- and B-correlations in 
the underlying irreducible structure, and, therefore 
they can be further reached by one Pauli loop. 

The two ^-independent vertex corrections V b and Vf 
are inter-related and can be computed by using the 
RFHNC summation technique. The series of one-body 
diagrams corresponding to V b is exemplified by Fig. |2j 
(a), whose sum is given by 



V b = V b {K) = V b (B) = 



h Vf 



(11) 



1 - h a V f 

where the Fourier transform h(k) of h(r) is defined as 



h{k) = p / dfe lk - r h(r) 



(12) 



and h = h(k = 0). 

The vertex correction Vf concerns the reducible 
fermionic diagrams, which in the boson limit [y = N) 
corresponds to unlinked bosonic diagrams (like diagram 
(la)), with a coefficient given by the next to leading order 
in N. 



A series of one-body diagrams contributing to Vf is 
displayed in Fig. |^ (b). The corresponding vertex cor- 
rection Si is given by 



Si = 1 - 2V b + 3V h 2 - 4K 3 



(1 + v b y 



(13) 



Besides Si , there is a second type of vertex correction, 
S2 , contributing to Vf . The diagrammatic representation 
of S2 is displayed in Fig. @ (c), and is, explicitly 



S 2 = 



dk 



h 2 {k)V f 



{2TTfpl-h 2 (k)Vf 



(14) 



The fermionic vertex correction Vf is made of vertices 
of the Si and S2 type, as schematically shown in Fig. || 
(d), and it is given by 



Si 



I + S1S2 (l + V b ) 2 + S 2 



(15) 



The vertex corrections V b and Vf correctly satisfies 
the normalization condition for the one-body distribu- 
tion function, <?i(l) = 1. The sum of the entire set of 
one-body diagrams is given by 



Si (1) = (1 + v b fv f + s 2 v f = 1 



Pair distribution function. 



(16) 



The calculation of g(r 12 ) can be done by using the 
FHNC summation technique for irreducible diagrams, 
renormalizing the vertices with the vertex corrections V b 
and Vf , depending on the nature of the vertex itself. The 
independent pair nature of the trial function reduces the 
set of allowed cluster diagrams, and, as a consequence, 
the structure of the integral equations results to be quite 
different from that encountered in standard FHNC p3[ . 
The pair distribution function is a sum of chain and com- 
posite diagrams. 

The lowest order diagrams contributing to g(r 12) are 
(i) the ideal Bose gas (uncorrelated) diagram, r„(l,2), 
and the two-body correlated diagrams, either (ii) lin- 
ear, or (iii) quadratic in h. They are displayed in Fig. 
|3|, together with their expressions, which include ver- 
tex corrections. They correspond to the set of dia- 
grams included in the lowest order constrained varia- 
tional (LOCV) method [^9| , when the Jastrow correlation 
is taken of the form f(r) = 1 + Vfh(r) and the vertex cor- 
rection (1 + V b ) 2 Vf, which gives the condensate fraction 
(see Eq. (p9|)), is set equal to 1. 

The lowest order bond, used in the chain diagrams, is 
represented by the solid line and corresponds to the cor- 
relation function h(r). After the inclusion of the vertex 
corrections, it gives the lowest order composite diagram 
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Vb = . x = . # + ■ ♦ # + • • • (a) 





FIG. 2: Diagrammatic equations for the vertex corrections. Equation (a) gives the bosonic vertex correction Vb- Equations (b) 
and (c) refer to the vertex corrections Si and 52 respectively. Equation (d) exemplifies the cluster series of the full fermionic 
vertex correction Vf. The solid lines represent the correlation function h(r), whereas the oriented lines indicate the exchange 
function l(kpr). 



B, K 



2(1 



V b ) 2 Vf h(r l2 ) 



X (r 12 )=2(l + V b ) 2 Vf h{r 12 ) 



(17) 



There are four types of chain diagrams, depending on 
the nature (B. K) of the external vertices. For symme- 
try reasons, N BK (1,2) = N KB (2,1) = N KB (1,2) and 
N BB (1, 2) = NxxiX, 2). The allowed chain diagrams are 
only those constructed with the lowest order bond and 
they are given by 



K 



B 



Vfh 2 (r 12 ) 



FIG. 3: Vertex corrected two-body diagrams which con- 
tribute to the distribution function g(r\2), together with the 
expression of the corresponding cluster terms. The vertex 
corrections to the top diagrammatical structure r u (l, 2) sum 
up to unity (see Eq. (|16[)) 



N BK (q) = N KB (q) = (1 + V h ) 2 V f 



N BB (q) = N KK (q) = {l + V b ) 2 Vf 



2 ' 



h 2 (q)Vf 
l-h 2 (q)Vf 

h 3 (q)Vf 
l-h 2 (q)V 2 



(18) 



The sum of all vertex corrected chain diagrams is 



N(q) = 2(1 + V b yVf(N BK (q) + N BB (q)) 



2(1 + V b ) 2 V f 



h 2 (q)Vf 
1 - h(q)V f • 



(19) 



The vertex corrected composite diagrams contributing 
to g(r\2) are given by 
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III. DENSITY MATRICES. 



X(r l2 ) = X Q (r 12 ) + Xi(n 2 ) + X 2 {r 12 ) , 



(20) 



where X (r) is given by Eq. ( |l7j ) and the Fourier Trans- 
forms of X\{r) and X 2 (r) are 



Xi{q) = 



X 2 (q) = 



h?Vf 



K*Vf 



I - h 2 Vf I - h 2 Vf ) q 



hV f 



hV f 



l-h 2 V 2 l-h 2 V? 

J J / q 



(21) 



where the convolution integral (\) q is defined by 



(f\s) q 



p+q 



.., . dppf{p) / dkkg{k) . (22) 
(2nypq J J lp _ ql 



The sum of all the chain and composite diagrams gives 
the exact expression of the two-body distribution func- 
tion g(r 12) in the thermodynamic limit, 



g(n 2 ) = 1 + N{r 12 ) + X(r 12 ) + 



(23) 



The leading term of the 0(1 /N) terms is independent 
of the coordinate r 12 , and has the form a/N, where the 
constant a is given by 



l-p / df 12 (N(r 12 )+X(r 12 )) 



(24) 



Notice that the constant a/N, contributes to S(0), and 
makes the normalization sum rule of Eq. (|j) to be cor- 
rectly satisfied, independently on the correlation function 
h{r) considered. The set of cluster diagrams giving a 
is constituted by the renormalized two-body Pauli loop 
(the exchange of the uncorrelated diagram r„(l,2) in 
Fig. ||) plus all those corresponding to its the convo- 
lution with 5(^12) — 1. Terms of the type f(r\ 2 )/N or 
higher order 1 /N terms can be easily calculated by using 
the proposed FHNC scheme, but they are not of interest 
for the present paper. 

The static structure function can be easily calculated 
from Eq. and it results to be 



S(k) 



1 + 2(1 



V b ) 2 V? 



h(k) 



1 - h(k)V f 



h 2 Vf 



h?Vf l-h 2 Vf t 



hV f 



hV f 



\l - h 2 Vf 1 - h 2 V 2 1 



(25) 



We compute in this section the one- and two-body 
density matrices of a bulk boson system described by the 
IPC trial function. They give us important information, 
such as the condensate fraction, the (paired) depletion 
momentum distribution and the pairing function. 

The one-body density matrix pi(l, 1') is defined as 



p(l,l')=p(rii')= (26) 
/ dr 2 dr 3 . . . dr N <S>* IPC (V \ 2,..., N)* IPC {1, 2,...,N) 



J dR\* IP c(R)\' 



and its related quantity, the one-particle momentum 
distribution is given by 



n(k) = (ala~) 



p j drive % "- r ^' p{r lv ) 
n Q 6{k)+n d (k) , (27) 



defining a condensate and boson-pair number distri- 
bution that can be independently evaluated in this IPC 
scheme. Notice that, with the above definitions, the mo- 
mentum distribution normalization is a sum rule on uq 
and rid(k) 



n 



dk 



(2tt) 3 p 



n d {k) 



(28) 



A. Condensate fraction and depletion contribution 

The term independent of r\y in the one-body density 
matrix, is the condensate fraction, and it is given by the 
product of the vertex corrections in the points 1 and 1' 



n Q = lim p(r X v) = (1 + V b ) 2 V f = 



(1 



(1 + V b ) 2 + S, 



(.29) 



The coherent depletion term nd(k) is given by a di- 
agrammatic series characterized by chains of bonds, al- 
ternatively of the B- and the if-type, starting with the 
-B-type from point 1 and reaching 1' with the -ftT-type. 
Each vertex has a Vf correction (the external points 1 
and 1' have a single vertex correction). The sum of this 
series is given by 



n d {k) 



h 2 (k)V? 



l-h 2 {k)V 2 

The normalization sum rule (Eq) is fully satisfied 



(30) 
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2' 



2' 



2' 



dk 



-n(k) = 



{l + V b f 



+ S 2 V f = 1 , (31) 



(2^)3 P v ' (1 + V b ) 2 + S 2 
and the depletion of the condensate is given by 

1 - n = S 2 Vf . 
B. ODLRO and pairing behaviour. 



(32) 



The study of the off-diagonal long range order and the 
calculation of the paring function Q are based upon the 
two-body density matrix, which is defined by 



p 2 (l,2;l',2') 



N(N - 1) 



(33) 



J dr 3 dr 4 . . . dr N V*(l', 2', 3, ... , JV)W(1, 2, 3, . . . , JV) 

The diagrammatical expansion of p 2 (l, 2; 1', 2') is char- 
acterized by the four external vertices 1, 2, 1' and 2'. The 
particular form of the IPC wave function excludes the 
possibility that more than one /i-bond reach any of these 
vertices. It follows that, in the thermodynamic limit, 
the three unlinked structures in Fig. |4| are the only ones 
which contribute to p 2 , with the result 



P2 (l,2;l',2')=p(l,r>(2,2') 

+n Q (n d {r 12 >) + n d {r V2 ) + P{r 12 ) + P{r V2 ,)) 



+P(ri 2 )P(rv2') + n d {r 12l )n d {r V2 ) + I- (34) 

where n d (r) is the reverse-Fourier transform of the de- 
pletion momentum distribution n d (k) of Eq. ([30j), and 
the function P(r), defined by 

n Q P(r) = X (r) + N BB (r) = X Q (r) + N KK (r) , (35) 
is the reverse-Fourier transform of the pairing function 

lonll 



1' 



1' 



1' 



FIG. 4: Diagrammatical structures contributing to the two- 
particle density matrix. The four vertices of each structure 
are renormalized by their vertex corrections, as in Fig. []. 
Each bond denote either no correlations or a generic linked 
structure. The fully uncorrelated diagram should be counted 
only once. 



• The first diagram on the left of Fig J corresponds 
to the term p(l, l')p(2, 2') in Eq. (p4|), which also 
includes the diagram with no bonds giving Uq. 

• Diagrams belonging to the other two structures, 
having one bond only, contribute to terms of the 
type non d and UqP. Those with two bonds give 
rise to the pairing and anti-pairing terms P x P 
and n d x n d respectively. 

• By using the HNC theory developed in section 2, 
one can easily verify that the leading 0(1/N) term 
of Eq. (|34|) is given by a/N, exactly as for the two- 
body distribution function. 

Simple inspection of Eq. ( j34|) shows that the diagonal 
part of p 2 (l, 2; 1', 2') coincides with the two-body distri- 
bution function g(r\ 2 ), namely 



p 2 (l,2;l,2) = ff (r 12 ) 



(37) 



The following independent-particle behaviour of the 
density matrix is found if the two particles 1 and 2 are 
far away from each other 



lim P2 (l,2;l',2')-p(l,l>(2,2'). (38) 

ri2,TV a /— MX3 

However, there is another interesting limit, which is 
associated with the ODLRO. This is when r and r' are 
far away from each other 



P W = 1 H l2nL2 = -VMfcXl + Mfc)) ■ (36) 
1 — h z {k)Vj 

The three diagrammatical structures of Fig. refer to 



direct terms of the cluster expansion of p 2 of Eq. (|33[) , 
namely those terms in which the coordinate 1 in ^ cor- 
responds to 1' in \E r *, and 2 to 2'. The corresponding 
fermionic diagrams have an exchange line joining 1 with 
1' and another one joining 2 with 2', which are not drawn 
in the figure, and that, in the bosonic limit, tend to unity. 
The terms, in which 1 and 2 are exchanged are of order 
l/N or higher. 

The various terms in Eq. ( |34| ) have been obtained as 
follows. 



lim p 2 (l,2;T ) 2 , ) = 

r ll'. I- 22'^°° 

(no + P(r 12 )){n + P(ry 2l )) 



(39) 



indicating that the pairing function P(r) is related to 
the superfluidity properties of the boson model system 



The two particle momentum distribution n(k\,k 2 ) = 
at at at ar ) can be computed by Fourier transforming 



the two-body density matrix |3lj, 32 1, namely 



n 2 (ki,k 2 ) 



(40) 



^ / dr 11 ,dr 22 ^r 12 e lfel -^' + ^^> 2 (l,2;l',2') . 
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Eq. 



garantees that the normalization property 



where 



dk 1 dk 2 f f r . 
——^n 2 {k u k 2 ) 



(1 



1, 

N' 



(41) 



is fully satisfied. 



By inserting the r.h.s of Eq. (34) into Eq. (|4C|), we get 
the result 



n 2 (k 1 ,k 2 ) = n{ki)n(k 2 ) + 

k 2 )P{k 1 )P{k 2 ) 
5{k x - k 2 )n d {ki)n d (k 2 ) + S(k 1 )S(k 2 )ap) . (42) 



" \ 1 — h Vf 



Notice that the distribution of the (k, —k) pairs results 
to be rid(k) 2 + P(k) 2 . This means that the distribution 
probability n d {k) provided by a wave function with no 
ODLRO is increased by a finite amount given by P 2 (k). 
A similar feature is obtained for the (k, k) pairs, whose 
distribution probability is increased by n 2 d (k). 



IV. ENERGY EXPECTATION VALUE AND 
EULER EQUATION. 

The energy expectation value is defined as 



E[h] 



(V) <T> 
N N 

The potential energy (V) /N is given by 



(43) 



N 



-p I <lrr(r)<i(r) 



(44) 



2 P 



dq 



(2tt) 3 p 



v(q)(S(q) - 1) 



2 P 



drv(r) 



which can be computed by using either the expression 
for the pair function g(r) give n in Eq. ( |23| ) or the static 
structure function of Eq. (P5|) . Both g(r) and S(q) de- 
pend upon /i(r), therefore (V)/N is a functional of h(r). 

The kinetic energy expectation value is defined as 



t(r) = -l-i^rjr)" 
2m r 

t(q) = e(q)h(q) , 



(47) 



and e((7) = h 2 q 2 /(2m) is the uncorrelated single parti- 
cle energy. 

The cluster expansion of (T)/N is similar to that of 
the pair distribution function. The two-body diagram, 
with only one correlation line (t(r\ 2 )) between the two 
external points 1 and 2, gives no contribution under re- 
integration. Therefore, one is left with the two-body cor- 
related diagram, having both t(ri 2 ) and h{r\ 2 ) as corre- 
lation lines, plus the composite diagrams given by t{r\ 2 ) 
dressed by Nll{tv2) chains, leading to the following ex- 
pression 



(T) 
N 



dq t(q)h(q)V 2 
{2n) 3 P l-h 2 (q)V j 



f 



dq 



(2tt) 3 p 



:(q)n(q) .(48) 



The above equation proves that the kinetic energy sum 
rule is exactly satisfied in the IPC model with the FHNC 
treatment. 

The total energy per particle E = (T) /N + (V) /N can 
be expressed in terms of the momentum distribution, 



E[n d ] = 



%-e(q)n d (q) + ^(0) 



(27r)V"' ~" W ' 2 

dq 



(2ir) 3 p 

1 f dq 

2 J (2nfp 
1 f _dq 



(n d (q) + P{q))v(<l) 



v{q) (n d \n d ) a 



(49) 



A. Applying the Ritz variational principle. 

The Euler equation is obtained by applying the Ritz 
variational principle to the energy per particle E[h]. 



TV 



h 2 (tt|vj|w) 



(45) 



The action of on (R]^) does not lead to three-body 
operators as in the Jastrow case, since 



+ t(ru)h(n m ) + ■■■ (46) 

(i#l)/(/<m) 



Sh(q) 



E[h] = 



(50) 



which gives rise to an integro-differential equation hav- 
ing h(q) as functional variable, derived in the Appendix. 

Since the total energy can be expressed in terms of the 
background momentum distribution n d (q) only, as shown 
in Eq. (|49|), then Eq. (|50|) can be written as 



5E 



dk- 



SE 5n d {k) 
Sn d (k) 8h(q) 



(51) 
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and from Eq. @ 



SE[n d ] = q 2 
Sn d (q) 2n 2 p 



+ n v(q) 1 - 



2n d (q) + 1 



2^/n d (q){l + n d (q)) / 
2n d (q) + 1 



^ 2^(^(1 + n d (g)) V J 



(2tt) 3 p 



[n d (k) + P(k) 



(52) 



and the variation of n d (q) with respect to h(k) is given 

by 

= 2V f F(k)S(q -k) + S^2h(k)F(k) , (53) 

where the explicit expressions of the functions f(fc) 
and SVf /5 h(q) can be found in the Appendix (se e eqs. 
( |A2| ) and ( |Al| ) respectively). Inserting Eq. ( |53| ) into 
Eq. (|5i"|), the Eq. (^ ) can be written as 



Going back to Eq. ( |5l| ) we see that one optimal solu- 
tion is A — 0. However, in Eq. (|59|) below we see that 
A may be non zero in order for the boson-pair number 
distribution n d (k) ~ 1/k so the compressibility sum rule 
can be obeyed. Hence we consider only this solution. 



B. The IPC Euler equation 

It is well known that the depletion component of the 
momentum distribution, n d (q) has a l/q singularity in 
the long wavelength limit |20f . The behavior n d (q) — > 
nomc/(2q) at small q, where c is the sound velocity, 
uniquely depends upon the existence of a condensate and 
of a finite value of the compressibility (2(], [2l[] . 

According to this important property, the long wave- 
length behavior of n d (q) can be written as 



nd(q) 



q 



d x q 



as 



9^0+, (58) 



and, consequently, h(q)Vf — > — 1 + q/{2d-\) as q — > 0. 
Using the above long wavelength behaviors, and devel- 
oping the r.h.s. of Eq. ([57|) around q — 0, one gets the 
solution 



SE[n d ] 



k 2 



-A[n d ] , 



8n d (k) 2t: 2 p 
where the ^-independent term AfnJ is given by 



(54) 



A[n d ] 



Vf{ax + C) 



I dkn d {k) (l + n d {k)) x , (55) 



with the constants a\ and C defined in eqs. (A3) and 
( |A2| ) respectively. Inserting the explicit expression of 
oE[n d ] 1 8n d (k) , given in Eq. ( |52] ) into eqs. ( j54| ) and (55), 
we get the following system of integro-differential equa- 
tions 



2A 



dk 



Vf(ax+C)J (2w) 3 p 
where 



n d (k)(l + n d (kj) = A ,(56) 



A[n d (q)] = e(q)+n v(q) ( 1 
+ (v\n d ) 



2n d (q) + 1 



2 y /n d (q)(l + n d (q)) i 
2n d {q) + 1 



,v\P 

2y/n d (q){l + n d (q)) V / , 



dk 



(2tt) 3 p 



v(k) (n d (k)+P(k) 



(57) 



A = -2 



dk 



(2tt) 3 p 



v(k)P(k) 



(59) 



By substituting Eq. (|59j) in Eq. one gets the fol- 
lowing IPC Euler equation 



e(q) + n v{q) ( 1 
+ (v\n) g 



2n d (q) + 1 



2y/n d [q)[l+n d {q)) > 
2n d (q) + 1 



dk 



(2tt) 3 p 



,vP) 

2y/n d (q)(l + n d (q)) V Jq 

v(k) (n d (k) - P(kj) = , 



(60) 



which is the a central result of the present paper. It is 
important to notice that the solution of Eq. (|60|), auto- 
matically satisfies Eq. (p6[) . Substituting the expressions 
of <zi and C into Eq. (p6|)7 and using the long wavelength 
behavior hoVf = —1, one gets 



dk 



(2tt) 3 p 



n d (k) =l + 2V f 



hoVf =1- V, ,(61) 
(1-Zio^/) 3 4 1 



which can be used to express the vertex correction Vt 
in terms of the condensate fraction no, 



4 1 



dk 



(2tt) 3 p 



n d (fc) 



4n n 



(62) 
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The same results, Vb = 1/2 and Vt — 4no follow from 
the vertex correction equations (|Tlf ) and (15|) and the 
depletion of the condensate given in Eq. f|32 ). 

Therefore, the IPC Eulcr equation (|5C) is consis- 
tent with the Gavoret-Nozieres singularity for the long- 
wavelength behavior of the momentum distribution of a 
quantum Bose fluid. 

The long wavelength behavior of the static structure 
function S(q), resulting from the solution of the Euler 
equation (BOh is given by 



eigenfunction, S(0 + ) = 0. Therefore, pair correlations 
alone are not sufficient to bring S^O" 1 ") from one to zero. 
However, we will show that for the IPC model, and its 
Bogoliubov limit reasonable potentials yield S(0 + ) ~ 0. 
Furthemore dS{q + )/dq | ~ l/(2mc), namely the long 
wavelength behavior of S{q) is very close to that required 
by sum rules pi|. 



C. Bogoliubov approximation 



S(q) -» S(0+)+q 



S(0 + ) 



dk 



(2tt) 3 p 



dS(0+) 
dq 

P 2 (k) > 



q^0+ , 



(63) 



Notice that for the ideal Bose gas S(q) = S(0 + ) = 1, 
whereas for the Jastrow ansatz, as well as for the exact 



It is interesting to understand under which approxi- 
mations our HNC theory recovers the Bogoliubov-type 
result. Since the sum of the chain diagrams corresponds 
to RPA theory, we expect that neglecting the composite 
diagrams X\ and X2 in the HNC summation lead to the 
Bogoliubov approximation |]l5| . Under this approxima- 
tion, and setting the condensation fraction no = 1, the 
total energy per particle of Eq. ( 49 ) is given by 0, [b3] 



E c hain[nd] — E un [rid] 

dq 



(27T)3p 



(64) 



{q)rid{q) + v(q) \n d (q) - \/ n d (q) +n d {qf 



where 



E un [rid] = ~v(0) 



(65) 



is the energy upperbound obtained with the uncorre- 
cted variational function ^(i?) = 1. 

Using Eq. (|64|) for the total energy, the Euler equation 
becomes 



e(q)+v(q) 1- 



2n d (q) + 1 
2y/n d {q)+ n 2 d (q) 



= 



(66) 



Ebo 9 = 2^(0) 



dq 



(n B og(q) (e(q)+v(q)) 



{2-Kfp 

v(q)\Jn B og(q) + n 2 Bog {q)) ■ (68) 



The static structure function S Bog {q) has a long wave- 
length behavior characterized by 5*(0 + ) = 1 — no and 
dS{q)/dq\ q=0 + = n /(2 v / E un ). 

The Bogoliubov approximation, and therefore our in- 
dependent pair model, is best for those interactions such 
that the scattering length a coincides with its Born ap- 
proximation 



which yields the Bogoliubov solution for the momen- 
tum distribution 



n Bog (q) 



Mq) 

uj 2 (q) = e 2 (q) + 2e(q)v(q) 



(e(q) + V(q) - u{q)) , 



(67) 



Anh 2 



drv(r) 



Airph' 



■v(0) 



(69) 



Notice that a = aeorn in the case of a <5-force inter- 
action of the type v(r) — An'h 2 ad(r)/m §). For such 
an interaction, the Bogoliubov approximation coincides 
with the Lee and Yang |b3 low-density expression 



n B og(q) goes as ^ E un / {2q) in the long wavelength 
limit. This Bogoliubov momentum distribution function Ely 
has been used as starting point in the numerical solution 
of the IPC Euler equation (|60|). Inserting n Bog (q) back 
into Eq. (|64|) one recovers the Bogoliubov estimate for 
the total energy per particle q. 



2ma 2 



Anpa 3 



(70) 



and the momentum distribution goes as q 4 at large 



12 



TABLE I: Gaussian potentials used in the calculations of 
the weakly interacting Bose system. The strength vq and the 
energy E un of the uncorrelated system (^(-R) = 1) are given 
in units h 2 /(2ma 2 ). The width a and the effective range r e ff 
are in units of the scattering length a. 

Potential no a a/aB r eff E un /p 

Gl 1.1778 x 10" a 11.256 095 -230.089 13.228 
G2 1.0028 x 10~ 2 5.6127 0.90 -51.369 13.963 
G3 9.1978 x 10~ 2 2.7887 0.80 -9.759 15.708 



TABLE II: Soft-sphere potentials used in the calculations 
of the weakly interacting Bose system and taken from Ref. 
|j7| . The strength vo and the energy E un of the uncorrelated 
system (^(R) = 1) are given in units h 2 /(2ma 2 ). The range 
parameter R and the effective range r e // are in units of the 
scattering length a. 

Potential «o R a/aB r e ff E un /p 

SS10 6.8167 x 1(T 3 10 0.880 -29.936 14.280 
SS5 6.3086 x 10~ 2 5 0.761 -4.9637 16.513 



V. RESULTS 

In this section we present and discuss the results ob- 
tained for the equation of state, the two-body distribu- 
tion function, the momentum distribution and the pair- 
ing function of a weakly interacting Bose system within 
the independent pair model, in correspondence with two 
different class of potentials. 

For low density systems, like atomic vapours, of in- 
terest in the BEC physics, the only quantity which is 
considered to characterize the interaction is the scatter- 
ing length a, namely its strength. For such systems, 
a <C r , where r = (3/(47r y o)) 1 / 3 is the average distance 
between the particles. Recently, the possibility of creat- 
ing Bose systems with relatively large values of pa 3 has 
been opened up J|5|, pq| . For pa 3 > 10~ 4 the shape of 
the interaction becomes important for quantities such as 
the energy per particle and the condensate fraction [j37| . 

In this paper we consider the properties of a Bose gas 
for pa 3 in the range 10 -6 — 10 _1 , or, equivalently for 
0.016 < (a/r ) < 0.75. The pseudo-potentials v(r) used 
are purely repulsive and have a range of interaction which 
goes from ~ O.lro to ~ lOro We have considered two 
different choices for v(r) 

(i) Repulsive gaussian (G) potential defined by 



v(r) = vq e 



(71) 



where a is the width of the gaussian. Table | gives 
the values of the parameters vq and a of three dif- 
ferent choices of gaussian potentials, labelled with 
Gl, G2 and G3. 

(ii) Soft-sphere (SS) potential defined by 



v(r) 



vq if r < R , 
otherwise 



(72) 



where we have taken R — 5a and R = 10a as in 
[[57J. The values of the parameters are reported in 
Table |fi} 

Tables | and [FJ give also the effective range parame- 
ter r e ff, the ratio a/aBom and the energy upperbound 



TABLE III: Energy per particle of the IPC model divided by 
the density, Eipc/p, at various densities, for the five poten- 
tials considered. The IPC upperbounds are compared with 
the low-density expansion estimates of Eq. (fro]). The ener- 
gies are given in units of Tl 2 / (2ma 2 ) and the density p in a -3 . 
In these units Elyo/p — 4-7T = 12.556 



p 




Ely/p 


Gl 


G2 


G3 


SSW 


SS5 


10" 


(i 


12.63 


12.63 


12.63 


12.64 


12.65 


12.68 


10" 


5 


12.76 


12.73 


12.76 


12.80 


12.78 


12.85 


10" 


4 


13.17 


12.92 


13.05 


13.22 


13.11 


13.31 


10" 


3 


14.48 


13.11 


13.48 


14.06 


13.64 


14.31 


10" 


2 


18.62 


13.20 


13.81 


15.01 


14.08 


15.56 


10" 


1 


31.70 


13.22 


13.93 


15.55 


14.24 


16.33 



E un provided by the uncorrelated (ideal) Bose gas. The 
scattering length a has been calculated by using a stan- 
dard procedures Q . Notice in the case of a delta-type 
potential, when a/a^ mn = 1 and Eg og = Ely, the up- 
perbound E un coincides with the lowest order term of 
the Lee and Yang expansion, Elyo = (fi 2 1 (2ma 2 ))47rpa 3 . 
Elyo has been proved to be an energy lower bound in 
the limit of pa 3 — > and under general conditions of the 
interaction, which are met by the potentials considered 
here @. 

The Euler equation ( |60|) has been solved numerically 
by means of a self-consistent iterative procedure. Few 
iterations were sufficient to reach convergence for all the 
cases considered. The optimal nd{q) has been found 
to have the n^mcj '(2q) behavior in the long wavelength 
limit. The results for the IPC energy per particle of 
Eq. ( |49| ) are displayed in Table [II in correspondence with 
the five potentials considered, and in Figs. || and |^ for 
the soft-sphere potentials. Ejpc can be compared with 
the upperbounds provided by E un of Eq. flS5[), given in 
Tables | and |l| and the estimates from the low density 
expansion, Ely, of Eq. (f70"l), which are independent on 
the shape of the potential. Fig s. M and || compare the 
IPC results with Esog of Eq. (p8|) and with the avail- 
able DMC calculations |37|]. The figures display also the 
DMC |0 and Jastrow [}40| results for the hard-sphere 
potential. 

The above results deserve the following comments 



(i) the IPC upperbounds are above the Bogoliubov es- 
timates for all the pseudopotentials and the den- 
sity values considered. Therefore, the terms of the 
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I i_l i i i_l i i i_l i i i_l i i i_l i i i_l 
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p 

FIG. 5: Energy per particle or equation of state of a weakly interacting Bose gas for the soft-sphere potential SS5. The 
IPC results (solid line) are compared with the low density estimates Ely (dashed line), those coming from the Bogoliubov 
approximation (dotted line), Eb os , and the DMC calculations of Ref. (crosses). The DMC (triangles) and Jastrow ^(J 
(squares) results obtained in correspondence with the hard-sphere potential are also reported for completeness. The energies 
are in units of h 2 / (2ma 2 ), and the density p in units of a~ 3 . 
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FIG. 6: Energy per particle or equation of state of a weakly interacting Bose gas for the soft-sphere potential SS10. The 
IPC results (solid line) are compared with the low density estimates Ely (dashed line), those coming from the Bogoliubov 
approximation (dotted line), Eboq, and the DMC calculations of Ref. |}7| (crosses). The DMC |^| (triangles) and Jastrow ^] 
(squares) results obtained in correspndence with the hard-sphere potential are also reported for completeness. The energies 
are in units of Tl 2 /(2ma 2 ), and the density p in units of a~ 3 . 



hamiltonian neglected in Bogoliubov theory give 
a repulsive contribution. Notice that, contrary to 
Eipc-, E>Bog is not an energy upperbound. The rel- 
ative differences between the IPC and Bogoliubov 
results, {Eipc — Esog)/ Epog, are reported on Ta- 
ble [f^. One can see that Eipc and Eb 03 get closer 
and closer when either a/aeorn — > 1 or for increas- 



ing values of pa 3 . In these two limits they both 
approach E un from below. 

(ii) the low density expression of Eq. ( f7(i| ) starts fail- 
ing already at pa? ~ 1CP 3 . At this density value 
and above it Ely is always larger than Ejpc- The 
next to leading order in the low density expansion, 
giving a term in p 2 log(p) pj"| , improves the total 
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TABLE IV: Comparison between the IPC model and the 
Bogoliubov theory. The table reports the relative difference 
E ' P e b Eb ° 3 x 10 3 for the five potentials considered. The den- 
sity is given in units of a" 3 . 



p 


Gl 


G2 


G3 


SSIO 


SS5 


irr B 


2.589 


12.34 


68.05 


19.11 


109.6 


1(T 5 


2.053 


10.80 


62.62 


17.03 


101.4 


10" 4 


1.143 


7.541 


49.54 


12.51 


81.32 


10~ 3 


0.355 


3.353 


28.19 


6.337 


47.82 


icr 2 


0.063 


0.836 


9.934 


2.451 


18.43 


10- 1 


0.010 


0.163 


2.713 


0.276 


9.003 



The arrows in Fig. |8J denote the bose wavevectors 
for the densities considered. For the typical values of 
the scattering length and mean interparticle spacing, 50 
A and 2000 A respectively, typical of BEC |Q, we find 
that e p w 1.1210" 11 eV and T b m 0.134^ K. 

To better show the q~ l behavior in the long wavelength 
limit, the figure plots qn c i(q) rather than nd(q)- Similar 
results have been obtained for the other four potentials. 
The extrapolated values at q — > ar e comp ared with 
the theoretical sum rule value dr~± = yJn^JC/A where the 
compressibility K, is given by 



energy estimates up to pa 3 ~ 10~ 4 (adding this 
correction term to Ely the low density expansion 
agrees with DMC |$7| within three digits). How- 
ever, it gives too much attraction at higher density 
values. For instance, in the units of Fig. ||, the cor- 



rection at pa 
is -11.37 



3 _ 



10~ s is -1.706, and at pa 3 = 10" 



dp \ dp 



We plot in Fig. the quantity 



mc 



(75) 



A = 



(76) 



(iii) there is an overall agreement of the IPC upper- 
bounds with the available DMC |37j results, even 
for the SS5 potential, which is the one with the 
shorter range of interaction. Contrary to Ely and 
to Eboqi both the IPC and DMC results are sand- 
wiched between E un and the lower bound Elyo in 
the full range of density values considered. 

(iv) the shape of the potential starts becoming impor- 
tant already at pa 3 ~ 10~ 4 . Therefore, in the re- 
gion of large a, the scattering length is not anymore 
sufficient to determine the interaction. 

The condensate fraction no is displayed in Fig. [?] for 
the potentials Gl, G3 and SS5 as a function of pa 3 , and 
compared with the low-density expansion (|, ||, |42) 



as a function of pa 3 for the potentials Gl, G3 and SS5. 
The compressibility sum rule is satisfied when A = 1. 
For the sake of clarity we also plot in the lower panel the 
following difference function 



5n d {q) 



nd(q) — — • 



(77) 



We have found that over a range of densities pa 3 < 
10~ 3 the compressibility sum rule is satisfied within three 
significant figures. 

The coefficient d-i can also be calculated directly, by 
using the following expression, which results from the 
terms in q^ 1 of the Euler equation (|60|). 



n 



1 




(73) 



One can see that the condensate fraction of the IPC 
model starts flattening at pa 3 w 10 -2 , as for the case of 
E/p. However, it remains below one up to the highest 
density value considered. 

The depletion momentum distribution n d (k) for the 
potential SS10 at three density values is displayed in Fig. 
^. The optimized distribution has interaction induced 
k (pair) occupation even at T = 0, reminiscent of 
the Fermi distribution. 

Following the Fermi analogue it is useful to formally 
define Bose wavevector q b and bose energy, e d and tem- 
perature Tf, through 



n d (qb) 



2 m 



(74) 



d\ = 



n oVo + $ j2^r p v{q)P(q) 



Vlv{q) 



(78) 



The coefficient d-\ obtained from the above equation 
and lim g ^o qn-d(q) coincide within three digits for all the 
cases considered in the present calculation. 

We have also computed the chemical potential defined 
by 



dE 

E + p— 
dp 



(79) 



The results obtained for the sound velocity and for the 
chemical potential are given in Table for the poten- 
tial Gl, which has the ratio a/aeorn closer to 1. The 
corresponding estimates resulting from the low density 
expansion of Eq. (|70|) , and given by 
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FIG. 7: Condensate fraction of a weakly interacting Bose gas as a function of the density for the potentials Gl, G3 and SS5. 
The low density expression of Eq. is also reported for comparison. 



TABLE V: Sound velocity c and chemical potential p at 
various density values for the potential Gl. The low-density 
results given in eqs. ( jil| ) are also reported for comparison. 
The density is givne in units of o -3 and the compressibility, 
mc 2 , and chemical potential, p, in units of h 2 /(2ma 2 ) 



p 


(mc 2 / ' p)ly 


(mc 2 /p)ipc 




{p/p)lPC 


10 _b 


25.36 


25.35 


25.28 


25.28 


1(T 5 


25.85 


25.67 


25.61 


25.52 


10" 4 


27.40 


26.11 


26.65 


25.92 


io~ 3 


32.31 


26.39 


29.92 


26.28 


10" 2 


47.82 


26.45 


40.26 


26.42 


10- 1 


96.88 


26.46 


72.96 


26.45 



Fig 111] are small, 0.0102, 0.0270 and 0.0551 for pa 3 = 
10~^10~ 4 , 10~ 3 respectively. The values extrapolated 
directly from S(q) and those computed from Eq. ( |63| ) 
coincide within three digits for all the cases considered. 

The pairing function P(r '12), as the one-body den- 
sity matrix p(rn'), is long-ranged. Its Fourier trans- 
form. P(q), has therefore a singular long wavelength be- 
havior, given by —nQmc/(2q), opposite to that of the 
momentum distribution n c i(q). Therefore, the function 
Fi(r) = P(r) + p(r) is a short-range function. This 
function is the key ingredient of the following important 
ODLRO property of the semi-diagonal two-body density 
matrix ElL W3] 



(mc 2 ) LY = (^L_) Snpa 3 + 16^- j , (80) 

(^>- 3 M^?)' (8i) 

are also reported in the table. 

The pair correlation function h(r) multiplied by the 
vertex correction Vf can be calculated from the momen- 
tum distribution rid(q), by using Eq. (^0|). Its Fourier 
transform is plotted versus q in Fig. |l0| at three different 
values of pa 3 for the potential SS10. 

Finally, the two-body distribution function g(r '12 ) and 
the pairing function P(r 12) are displayed in Figs. O and 
|l2| respectively, for the case of the potential SS10. Similar 
plots can be shown for the other potentials. The long 
wavelength limits of the static structure function, 5(0 + ), 
obtained in correspondence with the g{r\2) reported in 



lim p(l,2;l',2) - n (l + Fi(r 12 )) . (82) 

r' { — >oo 

One can easily verify that -Fi(O) = —1/2, as required 
by the Stringari sum rule on the two-body density ma- 
trix pl| . Such a behavior of F\{q) has important conse- 
quences on the overlap between the particle and density 
excited states |l4j . Moreover the sum rule relationg p to 
F\ [|l) is satisfied numerically to three significant figures 
|44j impliying gapless exitation. 

VI. CONCLUSIONS AND PERSPECTIVES 

In this paper we have used the variational theory to 
develop a number-conserving model for boson pairing 
based upon the Bogoliubov theory. This model consider 
an independent pair correlated wave function as trial 
function, which is characterized by a 100% condensation 
of (k, —k) pairs. 
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FIG. 8: Momentum distribution of the IPC model of a Bose gas interacting with the potential SS10 at several densities. The 
momentum distribution nd(q) is multiplied by q to make clear the q~ long wavelength behavior and divided by ^fp to fit the 
scale. The full circles in the lower panel indicate the Bogoliubov results, \fE un /2. 



The cluster expansions for the two-body distribution 
function g{r\2) and the one- and two-body density ma- 
trices from the IPC trial function have been derived, and 
a scheme of the HNC type has been developed to sum up 
the cluster terms to all order in a fully closed form. This 
novel HNC method maps the Bose system into a Fermi 
system, such that the fermions have as many flavours 
as the number of particles, and then it makes use of 
the Renormalizcd FHNC theory to compute the resulting 
cluster terms. 

As a consequence, the energy expectation value can 
be obtained with no approximation for any given pair 
correlation h(r), resulting in a true energy upperbound. 
Similarly, the normalization sum rules on g(r 12), p{ru>) 
and ^2(1, 2; 1', 2'), as well as the kinetic energy sum rule, 



are exactly satisfied. This cannot be achieved in Jas- 
trow theory, where truncations of the cluster series are 
intrinsically made. 

The Euler equation to get the optimal pair correlation 
h(r), obtained in accordance with the Ritz variational 
principle, has been derived. Since the total energy per 
particle can be written in terms of the momentum dis- 
tribution only, we could derive the integro-differential 
equation having the momentum distribution as variable, 
instead of h(r). We found that the Euler equation in 
nd(q) is not uniquely determined, and we selected the 
one consistent with a momentum distribution which has 
the Gavoret-Nozieres singularity in the long wavelength 
limit. We proved that such an Euler equation can be 
numerically solved by means of a rather simple iterative 
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FIG. 9: Compressibility sum rule. Upper plot: The quantity A of Eq. ([76]) is plotted as a function of the density for the 
potentials Gl, G3 and SS5. Lower plot: Difference function of the momentum distribution with respect to that given for the 
exact sum rule at low q a at different densites for the SS5 potential. 



self-consistent procedure, and that the resulting nd{q) 
has indeed the required long wavelength behavior. 

Calculations have been performed for five different 
pseudopotentials, three of the gaussian form, and two 
of the soft-sphere type, spanning a region of interaction 
range going from ~ O.lro to ~ 10ro and of strength a/ro 
from - 0.02 to - 0.8. 

We found that IPC results for the energy per parti- 
cle and the condensate fraction are in reasonably good 
agreement with the available DMC calculations, particu- 
larly for those potentials for which the scattering length 
a is closer to its value cJBom in Born approximation. In 
the limit of a/aBorn - ► 1 the low density expansion, Bo- 
goliubov theory, DMC calculations and IPC model all 
give very close results. However, even for o/osom ~ 0.7, 
IPC provides a quite realistic representation of the cor- 



responding Bose system. 

We were also able to compute other interesting quan- 
tities, like the sound velocity, the chemical potential and 
the pairing function. The last quantity is directly related 
to the off-diagonal long range order in the gas, and to 
its superfluid properties at zero temperature. Therefore 
the IPC model appears to be an excellent candidate to 
study boson pairing in a quantitative way. 

We also found that the range of validity of the low- 
density expansion is limited to pa 3 ~ 10 -4 . Accordingly, 
for values of pa 3 of the order or larger than 10 -3 , the 
results for the various quantities of interest do dipend 
significantly on the shape of the potential and not only 
on the scattering length a. 

The proposed theory can be easily extended to treat 
finite Bose systems, typical of trapped atomic vapours 
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FIG. 10: Pair correlation function of the IPC model. The table displays the Fourier transforms of the pair correlation functions 
h(q) vertex corrected by V/, obtained for the case of the potential SS10. 



produced to study BEC and macroscopic coherent phe- 
nomena. Work in this direction is in progress. 

The generalization to Fermi systems or to Fermi-Bosc 
mixtures is also possible, although less straightforward. 
In this case, contrary to the Bose case, one has to deal 
with the elementary diagrams. In the Fermi case, the 
generic point of a cluster diagram can be reached by up 
to four bonds (instead of two), two of the /i— type and 
two of the statistical type. This feature still garantees a 
closed set of FHNC integral equations, but with kernels 
given by four-body functions. 

A third, challenging perspective is the study of the ex- 
citation spectrum and of the superfluid properties. Such 
studies draw on the IPC physical picture of a number- 
conserving ground state with both a condensate and a 
coherent zero-momentum paired depletion; and on the 
HNC treatment developed in this paper, which allows 
for exact evaluation of the expectation value of a generic 
n-body operator. 



APPENDIX A 

In this Appendix we derive the Euler equation resulting 
from the Ritz variational principle expressed by Eq. (|5C|). 

Let us compute the functional derivatives of the vertex 
corrections S2 and Vt. The vertex correction VJ, can be 
expressed in terms of Vf. 



SS 2 
Sh(q) 

SVf 
Sh(q) 

with 



6h(q) 



2ho SVf \ 
(1 - hoVf ) 3 5h{q) J 



C 



F(q) = 



dk h 2 (k){l + h 2 {k)Vf) 
(2tt) 3 p {l-h*(k)Vf)* 

Hq)Vf 



,(A1) 



{l-h?{k)Vj) 2 
The solution of the above linear system is given by 



(A2) 
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SS 2 
Sh(q) 

Sh(q) 
a 1 



7r 2 p \ax + C 

-iL { 1 
tt 2 p \ai + C 

2h 1 



F(q) 



(1 - hoV f f V] 



f 



(A3) 



The functional variation of the expectation value of the 
kinetic energy results to be 
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+ 



6h(q) 
SV f 



TT 2 p ) 2m 



(A4) 



dk 



e~h\k)vf) 



Sh{q) \ V f mj (2n) 3 p (l - h 2 (k)V 2 ) 2 / 
which leads to 



5(T) 
6h(q) 



i r die 



k 2 h 2 (k) 



fp{l-h 2 {k)Vf) 2 



(A5) 



The functional variation of the potential energy expec- 
tation value (V) is given by 



( Vb(q) 



S{v) -^)^m +vM - vi )' (M) 



5h(q) 

where V a (q) and Vb(q) arc defined as 



Va(q) = v\ 



h 2 V 2 \ 



2 I (l-h V f ) 2 



+ (l + h 2 (q)V?) 



(A7) 
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and the constant term V± by 



r i f 1 + hpVf r dk h(k)v(k) F ^[2^ {q2 ~ 2Tl) + J^) +Va{q) ' Vl ) =0 - {A8) 

1 ~~ ai + C \(l-hoVf) 3 J {2n) 3 Pl~h(k)V f 

1 f dk h(k)v(k) Since cannot be equal to zero, then the Euler 

equation becomes 



(1 - h V f ) 2 J (2^) 3 p (1 - h(k)V f y 

(2irfp tfVf)^ l-h 2 Vf ) k —q 2 h(q) + [V a (q)-Xyh(q) + V b (q)=0, (A9) 

J (2lT) 3 p ( ' \ (1 - h 2 V?) 2 1 ! _ h 2 V 2 J J 

We can now write down explicitly the Euler equation J% 2 
resulting from Eq. (|0|) A = — T i + V i ■ ( A1 °) 
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The quantities V a (q), Vb(q) and A all depe nd on h(r ) 
in a highly non linear way. Therefore eqs.(A9) and ( AlC ) 
need to be solved by means of a self consistent iterative 



procedure. The equation ( |A9[ ) is equivalent to Eq. (|6 
as one can verify by expressing h(q) in terms of rid{q) 
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